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Abstract 

We present an algorithm for rigid body diffusion Monte Carlo with importance sam- 
pling, which is based on a rigorous short-time expansion of the Green's function for 
rotational motion in three dimensions. We show that this short-time approximation 
provides correct sampling of the angular degrees of freedom, and provides a general 
way to incorporate importance sampling for all degrees of freedom. The full impor- 
tance sampling algorithm significantly improves both calculational efficiency and 
accuracy of ground state properties, and allows rotational and bending excitations 
in molecular van der Waals clusters to be studied directly. 



PACS NUMBERS: 02.70.Tt, 02.70.Uu, 36.40.-c, 67.40.Yv 



1 Introduction 



Diffusion Monte Carlo (DMC) methods provide a powerful theoretical ap- 
proach to the study of weakly bound clusters of atoms and molecules. The 
favorable polynomial scaling of DMC enables clusters with large numbers of 
degrees of freedom to be studied, provided that suitable schemes for sampling 
the multi-dimensional wavefunctions exist. In the case of van der Waals clus- 
ters containing one or more molecules, there often exist several different time 
scales, deriving on the one hand from the monomer internal degrees of freedom, 
and on the other hand from the 'external' monomer translations and rotations 
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which are associated with the van der Waals modes between monomers. Typi- 
cally, there is a large time scale separation between the high frequency internal 
vibrations of the monomer constituents, and lower frequency intermonomer 
van der Waals modes. This time scale separation is often invoked to motivate 
approximations in which one effectively freezes the high frequency degrees of 
freedom when studying the role of the low frequency modes. When applied to 
intramolecular degrees of freedom using Diffusion Monte Carlo, this approxi- 
mation leads to an algorithm first proposed by Buch [1], commonly referred to 
as 'Rigid Body Diffusion Monte Carlo' (RBDMC), in which all monomers are 
treated as non- vibrating, rigid molecules which are free to rotate and translate 
in the presence of intermolecular interactions. 

The first implementation of RBDMC for molecular clusters by Buch employed 
a short-time Green's function for the free rotation of each molecular monomer 
in its instantaneous principal axis frame[l] . A convenient Gaussian short-time 

approximation for such monomer rotational motion was assumed in that work. 
A number of applications of RBDMC have since been carried out, all of which 
use a combination of this angular sampling for monomer rotation in local 
monomer principal axis frames, together with the conventional sampling of 
translational motions of monomers in the laboratory frame[l-16]. The rigid 
body rotations are carried out using either Euler angles, or quaternion alge- 
bra [16]. Considerable use has been made of this algorithm for water clusters, 
where the floppy intermolecular degrees of freedom are difficult to study with 
other methods for clusters larger than the water dimer [1,2,4,6,7,10,11,14-16]. 

Nearly all of the RBDMC applications to date employ DMC in its simplest 
form, namely unbiased Monte Carlo sampling in which the ground state wave- 
function is sampled directly without using any trial function. Exceptions to 
this include a calculation for HF-Ar„ for which importance sampling of transla- 
tional degrees of freedom was implemented [3], and a calculation for rotational 
excitations in clusters of SFgHenfl?]. The latter calculation represents the first 
application of the complete importance sampling algorithm for rotations and 
translations which is described in detail in this work. Such biased sampling 
proves useful since unbiased sampling, while adequate for ground states of 
relatively strongly bound or fairly rigid clusters, rapidly becomes less efficient 
for weakly bound systems. Importance sampling also allows for easier com- 
putation of quantities other than the energy (e.g. structural quantities) that 
would require descendant weighting[18,19] in unbiased DMC. Furthermore, it 
is also more amenable to study of excited states than are unbiased sampling 
techniques. 

Biased sampling in DMC is generally achieved with the importance sampling 
approach first proposed by Kalos and co-workers for translational motions [20]. 
The improved efficiency is dramatic for very weakly bound systems such as 
the clusters of helium and hydrogen, which are dominated by quantum effects. 
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Both pure helium clustcrs[21] and doped hchum clusters[22] have been studied 
with this approach. Excited states can be obtained from fixed node calcula- 
tions in which there is no explicit trial function but where nodal constraints 
are nevertheless imposed directly in an otherwise unbiased sampling proce- 
dure[23]. More efficient samphng of excited states within fixed node is often 
possible by using explicit trial functional forms. Exploration of nodal release, 
nodal optimization, spectral evolution and other approaches to go beyond the 
fixed node approximation is an active area of research[10, 24-34]. 

In this paper we first outline the standard rigid body diffusion Monte Carlo 
algorithm with unbiased sampling (Section 2). Our treatment includes a quan- 
tum mechanical derivation of the short-time Green's function for free rotation 
of a general asymmetric top (Section 2.1). For spherical tops, formal solutions 
for both the imaginary time, diffusive Green's function[35] and the real time 
Green's function[36] were established a long time ago. However for a general 
rigid body, i.e., an asymmetric top, formal solutions have been restricted to 
eigenfunction expansions [37,38], and no rigorous derivation of the expected 
Gaussian form for small angle diffusion has been presented. A formal solution 
for this is useful since it enables extension for rotation in the presence of exter- 
nal force fields (to which the quantum forces in biased diffusion Monte Carlo 
are formally equivalent) to be naturally and rigorously implemented. In Sec- 
tion 2.2 we show how the use of a rigid body diffusion Monte Carlo algorithm 
is motivated, with a detailed discussion of the factors influencing time step 
choice in a weakly bound system composed of both atoms and molecules. In 
Section 3 we then extend the rotational Green's function formalism to enable 
direct importance sampling of rotational degrees of freedom. Implementation 
of the full importance sampled RBDMC algorithm is described in Section 4. In 
Section 5 we provide some basic applications of this approach which illustrate 
its use. First, we demonstrate explicitly with detailed tests that such short- 
time sampling of the angular degrees of freedom will cover the entire angular 
configuration space of a rigid body (Section 5.1). Then we show with specific 
examples taken from studies of two weakly bound molecule-helium systems, 
i.e., SFg ^He„ and HCN ^He„, that full importance sampled RBDMC consid- 
erably improves the efficiency of calculations for ground states (Section 5.2), 
as well as enabling calculation of Van der Waals excitations in such systems 
(Section 5.3). 



2 Rigid Body Diffusion Monte Carlo with Unbiased Sampling 
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2.1 Theory 



The diffusion Monte Carlo method provides a stochastic approach to solution 
of the Schrodinger equation 

a'jm^Hm)). (1) 



This is accomplished by first transforming the equation to imaginary time, 



dr 



mir)): (2) 



and then solving the resulting relaxation equation with the use of a short time 
approximation to the Green's function for H. The formal solution to Eq. (2) 
is 

|*(r)) = e-^n*(0)). (3) 



Expanding the initial wave function in a complete set of eigenfunctions of 

\^m = i:cnm, (4) 



we obtain 

I^W)=E^-e"''"1V^n). (5) 



Choosing the energy scale such that = will guarantee that the ground 
state solution is achieved asymptotically: 



Mr)) l^o) (6) 



The solution to Eq. (3) can be expressed in position representation by making 
use of the short-time Green's function: 



^{Q', T + At) = 1 (g'le-^^^ig)^(g, T)dQ (7) 
^JciQ^ Q', Ar)^(Q, T)dQ. (8) 
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— * 

where Q refers to the set of coordinates specifying a single body, and Q to the 
entire set for the A^-body system. In an ah-atom calculation, these coordinates 
will consist of the 3 A^- dimensional cartesian coordinates of the A^-body system, 
i.e., Q = R. We shall use the notation R to refer to the three-dimensional 

— * 

(center of mass) cartesian coordinates of a single body, and R to the 3A^- 
dimensional cartesian coordinates of the A'"-body system. If m of these bodies 
are molecules and have to be treated as rigid bodies, we have Q = R,Q, 
where Q = ((/?, i}, x) denotes the Euler angles specifying the orientation of the 
principal axis frame (PAF) of a monomer in the laboratory reference frame, 

— * 

Q denotes the 3m-dimensional set of angular coordinates for m monomers, 

— * 

and R refers to the combined set of 3m cartesian coordinates of the monomer 
centers of mass, and the 3n atomic coordinates for the n — N — m atoms. 

Eq. (8) can be propagated in imaginary time r to find the asymptotic solution, 
Eq. (6), given an accurate short-time Green's function. The propagation is 
started from some initial guess for the wave function, '^{Q,t — 0), which 
is represented as an ensemble of multi-dimensional configuration points or 
"walkers" : 

*(g,T = o) = ^<5(g,-g). (9) 



Thus the algorithmical issue reduces to finding an accurate short time Green's 
function for a given Hamiltonian H , together with an efficient sampling pro- 
cedure for the multi-dimensional integral in Eq. (8). 

The simplest diffusion Monte Carlo procedure samples the wave function 
^((5,r) directly. This unbiased, or non-importance sampled DMC has been 
extensively used in the past for applications to both electronic and nuclear 
structure and energetics[18,23,39]. We now summarize the conventional im- 
plementation of RBDMC using unbiased DMC, for Hamiltonians describing 
the nuclear motion of N-body atomic or molecular systems, with 



H = f + V - Er. (10) 

Er is a reference energy whose significance will become apparent below. We 
assume for simplicity that the interaction potential is pair-wise additive. 



V^J^^i^iv^i^^j)^ (11) 

i<j 

where Vij is the vector between the centers of mass Ri and Rj of monomers i 
and j, respectively. 
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Most conventional solutions implicitly utilize the second-order, split operator 
decomposition [40] of the Green's function: 

G = e-^^' ^ e-^^-e-^^-e-^^- + 0((Ar)2). (12) 



For Hamiltonians involving only translational kinetic energy, i.e., 

f'^j:-D,Vl (13) 
1=1 



where Di = / 27711, the contribution e"^*^"^ to the short time Green's function 
is given by the solution to the diffusion equation. In the position representation 
this is 



s 3N/2 



G'(R^R',Ar)=^—) (14) 



n 



/ 1 \3/2 



g -(R'-Kf/W,AT^-AT[iViR',n')+V{R,n))/2-ER]_ 



Note that this corresponds formally to rigid bodies moving with fixed relative 
orientations specified by the set of Euler angles {il^}. In general, this would 
not be a physically useful representation. However for a single molecule moving 
in an atomic cluster, e.g., in a rare gas cluster, this can sometimes provide a 
good representation, when the molecule is sufficiently heavy that the molecular 
rotational kinetic energy is small compared to the rare gas kinetic energy 
contributions [41]. The potential term in Eq. (12) is diagonal in the position 
representation: 

{Q\e-^^^\Q') = e ^""^{QIQ'). (15) 



In the position representation, Eq. (12) then becomes a product of non-local 
single particle diffusive terms from the kinetic energy, and a local term from 
the potential energy. Depending on its sign, the latter can lead to exponential 
growth or decay of the amplitude ^{Q) in Eq. (8). 

For Hamiltonians involving only rotational kinetic energy, we have for a general 
rigid body, i.e., an asymmetric top, 

p,a 



where dpa = 1/2/po,. Ipa is the principal axis moment of inertia for parti- 
cle p relative to its principal axis a — x,y, z, and Lpa is the corresponding 
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component of the angular momentum operator in the pth PAF: 

f h d 

Lpa^ x,y,z. (17) 

I d(ppa 

Here (p = {(px, <Py, <Pz) specifies the angles of one-dimensional rotations around 
the molecular principal axes x, y, and z, respectively. The rotational Green's 
function for a general rigid body should formally be expressed as G'^{Q, — > 
n', At). However, applying the second-order short-time decomposition 

P 



where 

g-T?AT ^ ^-d^Ll^Ar ^-dytlyAr ^-d.Ll^Ar _^ 0{At'^), (19) 



leads immediately to the following short-time factorization into one-dimensional 
rotational Green's functions: 

G{$ ^ 0', At) = n ^ ^'^ ' (20) 
p 



< 0;|e-'^-^-^n0, > . (21) 

Each of the one-dimensional rotational Green's functions can be evaluated by 
inserting a complete set of one-dimensional angular momentum eigenstates of 
the angular momentum along the corresponding molecular PAF. This is given 

by 

El^ra)(^ral = l' (22) 
m 



where 



LZ^) = e+^"^'^^«, (23) 



with eigenvalues 

4„|L- ) = m^|L- ), m = 0, ±1, ±2, ... (24) 
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Then, e.g., for o; = x, we readily obtain 



27r m 



where we have used the scaled constant dpa = h^dpa = /2Ipa for particle p 
and principal axis a. We shall see below that in the short time limit, Ar — > 0, 
this constant is identical to the rotational diffusion constant for this degree of 
freedom. Note that it is also equivalent to the molecular rotation constant Bp^. 
As discussed in Section 2.2 below, the magnitude of this relative to the masses 
of the n = N — m atoms, of the m molecules, together with the strength of 
the interaction potentials, will determine the time step for a given imaginary 
time propagation. 

It remains to be shown that Eq. (25) reduces to the expected Gaussian form 
in the small Ar limit. We do this by showing the reverse, i.e., by expanding 
the Gaussian form 



in eigenfunctions of the one- dimensional angular momentum operator L ~ 
—ihd/d(f): 

±oo -1 

9i<t>) = E ^--7|^^^"'^'- (27) 



n=0 



Hence 



i2 



For Ar << 1, we can extend the limits of this integral to ±oo, allowing 
evaluation to yield 



= ^=e-" "^^^ (29) 



Therefore we make the identification 



±00 -1 A 1 

^27r V47rdAr ' ^ ^ 



n=0 
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Replacing in Eq. (30) by {(ppx — (pprj) leads to the desired Gaussian form of 
the sum in Eq. (25), i.e., 

27rtr V^TTdAr ' ^ ' 



This yields the following short-time factorized limiting form for G\: 

At^O 1 -p-j- 1 _UA, ^2 

(i^^A^y (d^ 



The constants dp^ are now seen to act as one-dimensional diffusion constants 
for each of the rotational degrees of freedom in the p-th molecular PAF, as 
noted above. 

We note that Eq. (32) is formally overcomplete in angle space, giving rise to 
a total integrated angular volume of (Stt)^ instead of the requisite 47r^. This 
is a consequence of the short-time factorization into three one-dimensional 
rotations, Eq. (19). For a spherical top this can be avoided when the Eulerian 
short-time limit described below is employed. A comparison of the two short- 
time propagators for a spherical top is provided in Section 5, where we show 
that the formal violation of the angular configuration volume deriving from 
the product of one-dimensional rotations does not cause problems if the time 
step Ar is sufficiently small. The short-time factorization can be improved 
upon if necessary by making use of the split operator decomposition[40], but 
we have not found this necessary in any applications to date. 

For the special case of a spherical top {dpa = dp,a = x, y, z), the exact form of 
the Green' s function is well known[36], as is its short-time diffusive limit [35] 

/ . \ 3/2 



where ( is the angle of rotation about an arbitrary axis. This short-time dif- 
fusive limit was already derived by Furry in 1957[35]. It follows directly from 
taking the small small Ar limit to the imaginary time transcription of the 
infinitesimal propagator for a spherical top which was derived rigorously by 
Schulman (see Eq. (4.9) in Ref.[36], with F = C, and dp = h^/2I). Eq. (33) cor- 
responds to the Eulerian description of rotation of a rigid body as a rotation 
through a single angle ( about an arbitrary axis n(^[42]. It can be implemented 
in a Monte Carlo sampling procedure by first choosing a randomly oriented 
axis (uniform sampling of azimuthal angle (p on the interval [0, 2tt) and of 
cos(^) on the interval [—1, 0]), and then sampling the angle of rotation ( about 
this axis from a Gaussian distribution of width J3 x 2dpAr. The origin of the 
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factor of 3 is described in Appendix A. Related Eulerian rotational sampling 
has been employed in classical Monte Carlo rigid body simulations, but with 
the Gaussian sampling of the rotation angle ^ replaced by uniform sampling 
inside a fixed interval [—A(q/2, ACo/2][43]. This satisfies the conditions of re- 
versibility, detailed balance, and accessibility of all orientations, but does not 
correspond to the true Green's function, Eq. (33). It will therefore not neces- 
sarily yield the correct energies and expectation values in a quantum Monte 
Carlo calculation. 

When the Hamiltonian involves both translational and rotational kinetic en- 
ergy terms, the Green's functions G* and can be combined: 

G{R, n R', n', At) ^ G'{R R', Ar)G'"((^ ^ (f', At). (34) 

A multi-dimensional configuration (commonly referred to as a 'walker'), now 
carries both translational and angular coordinates. Sampling of the rotational 
and translational degrees of freedom are carried out by independent diffusive 
moves in these coordinates, as described by the factorization in Eq.(34). How- 
ever, there will be coupling of the translations with rotations via the potential 
energy term V{R,fl) in Eq. (12). Details of the stochastic evaluation of Eq. 
(15) is described in Section 4. Use of Eq. (34) is usually referred to as "rigid 
body diffusion Monte Carlo" (often abbreviated by its acronym RBDMC), 
because of the implicit representation of the full A^-body kinetic energy by a 
sum of rigid body rotational contributions together with the sum of the center 
of mass translational terms. 

2.2 Time step considerations 

The time step parameter used during the Monte Carlo simulation is deter- 
mined by two factors. First, as pointed out by Suhm and Watts [18], when all 
other factors are equal, the time step is controlled by the smallest mass. In- 
deed, for light particles, the diffusion coefficient is large and the time step has 
to be small enough so that the particle does not move too far in a single step. 
But this is not the only factor. The effect of the potential can also be impor- 
tant. For example, strongly bound particles must be sampled with a smaller 
displacement, and hence with a smaller time step, to avoid the majority of 
Monte Carlo configurational moves ending in destruction of the configuration. 
In the general case of a mixture of light and heavy particles, the time step has 
therefore to be chosen not only according to the smallest mass value, but also 
taking the effect of the potential into account. The presence of heavy, strongly 
bound particles may require the use of a shorter time step. For example, in 
the case of the water dimer, Gregory and Clary [44] use time steps of 0.5 - 
5.0 au for the all-atom calculations. Those small values are needed because of 
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the presence of high frequency internal modes of the monomers. When using 
RBDMC, they were able to increase the time step up to 30 - 60 au, because the 
introduction of the rigid body approximation eliminated these strongly con- 
fined degrees of freedom. This decrease in computational time step is generally 
cited as the primary motivation for an RBDMC calculation. 

The situation with regard to time steps is somewhat more complicated for 
helium clusters, where a wide range of time scales is operative. For pure he- 
lium clusters, the time steps used vary typically from 1000 to 10000 au [41,45]. 
Introduction of a non-rotating, rigid SFe impurity requires a time step reduc- 
tion to 250 - 500 au [41,46] even though the mass of SFg is much bigger than 
that of He. This reduction of time step is due to the fact that the helium- 
SFg interaction potential is considerably stronger than the He-He potential: 
the larger time step will lead to inefficient sampling of the He-SFg relative 
motion. When the molecular rotation is included, the situation becomes quite 
complicated as one now has the additional intrinsic time scale of the molecular 
rotation. Thus one needs to compare the corresponding diffusion coefficients 
dp and Dp for each species p included in the simulation, and to assess both the 
relative roles of the masses/moments of inertia (measured via d and D), and 
the confining role of the interactions in all degrees of freedom. For the systems 
studied in this paper, namely SFe and HCN in clusters of ^He, the defining 
parameters are given in Table 1. The time steps required to eliminate any 
time step bias for these two systems arc of order 20 - 40 au for SFg, 10 - 30 au 
for HCN) in implementations where attempted moves involve all particles be- 
ing moved together. Somewhat larger time steps can be used when attempted 
moves involve only single particles (e.g., ~ 50 au for SFe in ^He„). The latter 
is preferred for all except very diffuse systems. For molecules interacting with 
helium, in general all-atom potentials are not available, and this is a second 
reason why RBDMC is essential here. 



Table 1 

Diffusion coefficients and potential characteristics for SFg and HCN in ^He clusters. 
For SFg, V*"'' corresponds to an adiabatic barrier for a helium atom to move between 
sites. It corresponds to the minimum of the potential along the C2 symmetry axis. 
For HCN, V*"'^ is defined as the maximum of minrV {r,6) for 6 varying from to 
TT. It corresponds to the minimum of the potential for 9 = it. 



Species (X) 


D { 


au) 


d (au) 


^He-X m ^ 


Vf tx in K 


Eo in K 


He 


6.80 


10-^ 




-10.95 






SFg 


1.86 


10-^ 


4.15 lO"'^ 


-83.86 


-63.13 


-37.4 


HCN 


1.01 


10-5 


6.73 10-6 


-42.40 


-30.14 


-13.9 
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3 Rigid Body Diffusion Monte Carlo with Importance sampling 



We consider first the conventional form of importance sampling in cartesian 
coordinate systems. This starts with the introduction of a distribution function 
T{R,t) which is biased by a trial function ^t(-R): 

J^{R,t) = ^{R,t)'^t{R). (35) 



Taking its imaginary time derivative and using Eqs. (2) and (10) with T = T* 
(Eq. (13)), leads to 



^^^lll = [DV^ -DV -F-DF-V - El{R) + Er]T{R, t) 

^-H*T{R,t) (36) 

where El{R) is the "local energy", given by 

El{R) = ^T{Ry^H^T{R) (37) 

and F{R) is referred to as the "quantum force", given 

F(^) = 2Vln*r(^) (38) 



For simplicity, we have set Dp — D for all p — 1...N. The solution to Eq. (36) 
is then 



J^{R', At) = J G\R R', At)J^{R, 0)dR, (39) 

where the translational importance sampling Green's function & is formally 
defined by 

= ^T\R)G'iR R', At)^t{R')- (40) 



The short-time solution for G* is well known [39]. It is usually derived within 
a Fokker-Planck formulation [4 7]. We derive here an alternative solution using 
operators which provides a introduction to the subsequent new derivation of 
the full translational and rotational importance sampled Green's function, G. 
& is initially cast in the position representation, 

G\R^R',At) = (i?'|e-^"^'|^). (41) 
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The local energy contributions can then be immediately evaluated: 



G\R R', At) = e-^^[(^i(-^')+^^L(«))/2-EK] 

^^,|gAr[W^-DV.F-DF.V]|^^_ (42) 

We now replace by the operator [P — —iTiV) and use the identity 

F -V -^^P ■ F -V ■ F (43) 
n 

(see Appendix B for proof), to obtain 

G\R^R',At) = e-^-[(^^(^^')+^^(^))/2-^H](^'|e^^[-H^^^'-i^^-^]|^).(44) 



Inserting a complete set of momentum states with p representing the 

eigenvalues of the momentum operator P, allows the last factor in Eq. (44) to 
be rewritten as 

where we have used 

= (27r^)-'/'ei^^-^. (46) 



Substituting k — p/h, and evaluating Eq. (46) using the standard integral 

j d^k 6-'^^'+**^^ = (^^) e-^'/^", (47) 

leads to the well known result, 

^-AT[{Ei^{R')+EUR))/2-ET]^-[R'-R-DF{R)AT]y4DAT_ (^^^^ 



Now we expand the discussion to include the rotational kinetic energy T^, and 
to incorporate an orientation dependent interaction potential, V{R,fl). The 
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biased distribution function becomes 



(49) 



Identifying the quantum forces for one- dimensional angular motion around 
each of the principal axes, a — x,y, z, for each monomer p, 



d 



pa 



d<t>, 



(50) 



pa 



and the angular momentum operators Lp in the p-th PAF, Eq. (17), leads to 
the generalization of Eq. (35) for a general rigid body Hamiltonian Eq. (10), 
with f = T* + f"" and interaction potential V{R, Q): 



dJ^{R,n,T) ^ 2 



dr 



■ [DV^ - DV ■ F - DF - V 

\ d d 

~l~ / ^ „ ,„ / ^ \ ^paTZl fia diafpa , 

9(PL ^ [ ocppa dep. 



d(t)p 



pa ^Ypa pa V. ^H^pa 

-El{R, n)+ER]J^{R, Qr) 
-HT{R, Qr). 



pa 



(51) 



This can be simplified using Eq. (43) and its analog for angular motions (Ap- 
pendix B), 



f if f _9fa 

J a r% . J- ^aja 



d(t)a 



(52) 



to arrive at 



dJ^{R, Cl, t) 
97 



D 



i£p.p-J2dpaL 

pa 



2 

pa 



(53) 



-ih ^ dpaLpafpa - El{R) Er 



pa 



J^{R,fl,T) 



^-Hj^{R,n,T). 



(54) 



The short-time Green's function for H is put into position representation using 
Eqs. (34) and (20), and then making use of Eqs. (12) and (15): 



G{R, n R', n', At) 
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^-{R'-R-DAtF{R))'^/4DAt ^-AT[{V{R',n')+V{R,n))/2-Eii] 

n (0^^|e-['^^«'^^«+^'*'^^«'^^«-'^^«]^^|0p«) (55) 

pa 

Inserting a complete set of one-dimensional angular momentum cigcnstates, 
Eq. (23), in each of the terms inside the product on the right hand side as 
before (Section 2) leads to, e.g., for a = x, 



with dpa = h^dpa as before. This can be shown to reduce to a Gaussian 
form in the small At limit, via the same procedure as employed in Section 2. 
Thus, replacing in Eq. (30) now by the angular displacement including 
the contribution from the angular quantum force, e.g., for a — x, one has 
i'Ppx ~ (f>'px ~ dpxArfpx), leads to the desired Gaussian form of the sum in 
Eq. (56), i.e.. 



^—dpxm?AT+im{<l>px—4>'px~'ipx 



L ^-dpxm'^AT+im{<ppx-4>'px-dpxATfix) ^^3*° 

^ g - {4>px-4>'px -dpx ^rfpxf/idpx At ^ 

^jA'Kdp^Ar 



We thereby arrive at the full, short-time factorized, limiting form for the trans- 
lational and rotational importance sampled Green's function for a set of in- 
teracting rigid bodies: 



G{R, Q ^ M', At) ~ L_^e-^R'-R-D^rm)?/^DAr 

g-AT[(l/(J?'4')+V'(R,f1))/2-i?fl] ^ 
1 

g-[(</>p„-</-;,„-dpaAT/pa)V4<ipaAT] ^gg^ 

■pa (47rAr(ipQ,)^^^ 

This may be generalized to multiple values of by replacing the first factor 
with 

/ 1 X 3Af/2 / 1 \ 3/2 

^ ^ j ^-{R:^-R^-DpATFp{Rp)Yl^DpAT _ (^gg^ 



Details of the implementation of Eq. (58) are presented in the next Section. 

We note that, for the specific case of diatomic molecules, Lewerenz has pro- 
posed another method to implement rigid rotor motion, namely the method 
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of adiabatic constraints [48]. This method is attractively simple and easy to 
implement for ground state studies of linear molecules [48-50], although it 
becomes considerably more involved when applied to more general rigid bod- 
ies [18]. Importance sampling of rotational degrees of freedom would be prob- 
lematic within such an approach however, since it is not clear how to combine 
this with the adiabatic constraints. In particular, it is not obvious at what 
stage the quantum forces should be introduced, i.e., before or after imposition 
of the adiabatic constraint. Recently, Sarsa et al have shown that classical 
constraint dynamics can be used to generally impose bond length and bond 
angle constraints, in both biased and unbiased DMC[34]. 



4 Mixed frame and fixed frame implementations 

In order to demonstrate Monte Carlo solution of Eq. (8) with Eqs. (34) 
and (58), we consider a system composed of m rigid rotors and n atoms 
(n + m = A^). The Monte Carlo propagation over a step Ar requires then 
3(m + n) translational diffusion moves, 3m rotational diffusion moves, and the 
local exponential growth/decay operations. The translational Green's func- 
tion, Eq. (15), consists of a product of a diffusive term deriving from T*, the 
translational kinetic energy, and a branching term deriving from the potential 
energy, V . Two kinds of moves are thereby made on an initially established 
ensemble of mult i- dimensional configuration points or "walkers", Eq. (9). The 
first kind of move is a diffusion from Rj to new positions i?^, implemented by 
sampling the vector Rj — Rj from a 3iV-dimensional Gaussian distribution. 



of width y 2DpAT in each of the 3A^ cartesian coordinates. The second move 
is a "branching" which modifies the weight of each walker in the ensemble 
by the exponential factor e""^"^'^^^^'^"*"^^^^)/^""^^!. This can be implemented by 
i) replicating/destroying configurations, or by ii) maintaining and updating 
a continuous weighting factor associated with each configuration, or by iii) 
a combination of weights and replication [47,49]. The role of the reference en- 
ergy Eji now becomes clear, i.e., it provides a means to reduce the fiuctuations 
deriving from the branching term. 

Propagation with G*, Eq.(48), is performed by a simple modification of the 
unbiased propagation described above, namely that the diffusive moves are 
modified to include the contribution from the quantum force F{R)- The co- 
ordinates of each particle are therefore updated according to, e.g., for Rp^, 




'px 



( 



dMR) 

dxp 



/MR) 



(60) 



16 



where Ap^; is a random number sampled from a Gaussian of width y^2DpAr. 

The associated modification of the branching term consists in using the local 
energy instead of the potential, as shown in the exponential factor 
q-At{{El{R')+El{R))/2-Er] ^YiQ update of the weights and/or the replication 
of walkers. 

The implementation for the rotational Green's function, unbiased G" and bi- 
ased G''\ is similar. In the unbiased formalism, diffusive rotational moves con- 
sist in sampling the three rotation angles for each rotor from a 3m-dimcnsional 



Gaussian distribution of width ^2dpa/^T. In the biased version, these rotation 
angles also include the quantum force : 

= A,. + 2d,.Ar [^^^1^] /MR, (61) 



For the full Green's function, both translational and rotational diffusion moves 
are made. We now give a more detailed explanation of the various possible 
implementations of this in different frames of reference. 



4.1 Mixed Frame 



In the general case, i.e. for more than one rotor, it is necessary to use differ- 
ent frames for the rotation and the translation moves. The rotation must be 
performed in the principal axes frame, PAF of each rotor, the translations are 
done in the laboratory frame, SF. The evaluation of the potential and of the 
trial function have to be done in consistent frames. At each time step in the 
random walk, three rotations around the principal axes of the molecule are 
made. During the time step r — Ar — > r, the PAF (PAFt-_at) is thus rotated 
into a new PAF: (PAF^). The rotation matrix describing the attempted move 
is : 

^M=4(0D-4(</'p-4(C) (62) 

where Ra{4>) represents a rotation of around the a-axis, and 0^ denotes the 
angle moved around this axis in time r. The corresponding relation between 
the coordinates in the new PAF and in the old can be derived using the 
following identity: 

= 4.(01) •^,'(0;)-^x(0D- (63) 

Here y {z") is the new orientation of the y [z] axis after the first (the first and 
second) rotation(s). This identity is easily verified by transforming each of the 
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rotations appropriately [51], e.g. 

R^>{<t^l)=k,m-Ry{<t^l)-Rlm 



(64) 



The matrix representation of Eq.(63) is given by: 
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— sin (f)z cos (f)z 


1 


v) 


PAFr 


^ 
^10 


\ 


sin 


(j)y COS 4>y ^ 
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(65) 



(66) 



If the initial orientation of the PAF is given by 3 Euler angles ((/9o, '&o, Xo) with 
(fo e [0,7r), t^o £ [0, tt] and xo G [OjTt), Eq. (66) can be used to compute the 
relation between the coordinates in the PAF and in the SF. Then 



i=l 



PAFr 



y 



V{t) 



SF 



with 



SF 



(67) 



7^" = 



cos </7o COS ■j?o COS Xo — sin (/7o sin xo 
— cos </7o cos 'j?o sin xo ~ sin (/7o cos Xo 
cos(/?o sini?o 

sin (/7o cos 'i?o cos xo + cos (/7o sin xo — sin 'Oq cos xo 

— sin (/7o cos ■j?o sin xo + cos (/?o cos xo sin 'j?o sin xo 

sin (/7o sin i?o cos ■j?o 



(68) 
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In an unbiased "pure" DMC calculation, (p^, <i>y and 0^ are sampled from a 
Gaussian distribution which is centered at zero and has width v^2rf^^Ar. In an 
importance sampled calculation, we have to add the quantum force, according 
to Eq. (61). Rotation around x modifies the orientation of the y principal axis 
of the rotor. The second rotation is thus made around the new orientation of 
this axis, followed by the third rotation around the subsequent new orientation 
of the z principal axis. Due to the small values of these angles, we neglect the 
non-commutation of these 3 rotations. 

In the Mixed Frame formulation, the SF cartesian coordinates for the centers 
of mass of all particles arc used at every time step in the scheme. The use 
of the V matrix allows us to compute the coordinates of the relative vectors 
between the atoms and the rotors, in the PAF of each rotor. These relative 
vectors are needed for evaluation of the potential term and of the trial wave 
function. 

As noted above, the Mixed Frame formulation is needed whenever the system 
contains more than one rotor. It is also required whenever the Euler angles 
defining the orientation of the PAF in the SF are needed. These angles are 
defined by the V matrix. Knowledge of these angles is necessary for calculation 
of excited states within a fixed node approximation whenever the nodes are 
imposed in a SF frame, because the trial function is then an explicit function 
of these angles. In the particular case of a trial function that depends only 
on the second Euler angle , e.g. = 1 1,0,0) oc cos'd — V33, computation 
and storage of only the last column of V is sufficient. Section 5.3 provides an 
illustrative example of calculation of internal bending excitations for He-HCN 
using this Mixed Frame formulation. 



4-2 Fixed Molecular Frame 

When the system under consideration contains only one rotor, the Monte Carlo 
propagation can be performed in the PAF of this rotor. Like in the previous 
section, at each time step in the random walk, three rotations around the 
principal axes of the molecule are made. The corresponding relations between 
the coordinates are given by equation (66). 

In order to perform the translational moves, one continuously follows the 
molecular orientation. Then at each time step (r), the atoms and the center of 
mass of the rotor are defined by their coordinates in the molecular PAF.^- The 
matrix TZ^'^'' is used to compute the coordinates in the new PAF using their 
coordinates in the former one (PAF,-_At)- The translational moves are then 
made according to Eq. (60), using the coordinates in the old PAF. Vectors 
describing the position of the atoms in the two PAF's are used to compute 
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the atom-rotor interaction potential and the trial wave function, before and 
after the step. 

Calculations for SFg-He system using the implementation described in this 
section leads to the same results as those obtained with the Fixed Laboratory 
Frame. One disadvantage of this Fixed Molecular Frame scheme is that it does 
not allow for easy computation of the Euler angles defining the orientation of 
the PAF with respect to the SF frame. 

4-3 Fixed Laboratory Frame 

In the case of a spherical top molecule, Eq. (17) is applicable in the space-fixed 
laboratory reference frame (SF) as well as in the PAF, because any choice of 
orthogonal frame diagonalizes the inertia tensor. The rotational moves can 
therefore be made around the axes of the SF frame {X, Y, Z) . Although the 
rotation moves can be made without knowing the orientation of the rotor at a 
given time, it is necessary to keep track of the molecular orientation in the SF 
frame because computation of the interaction with the other rotors and atoms 
(Eq. (11)) requires knowledge of the orientation of the preassigned reference 
directions of the spherical top (e.g. the S — F bonds in the case of SFe He„ 
clusters) . 

Rotational moves, performed according to Eq. (63), are combined with trans- 
lational moves of the n atoms and of the m center of masses of the rotors, 
made in the cartesian coordinates of this Laboratory Frame (SF), within the 
conventional approach (Section 4 above). This fixed Laboratory Frame scheme 
has recently been implemented in a study of rotational excitations of SFg in 
hehum clusters [17], and is also used in the SFe applications presented here. 



5 Applications 

5.1 Angular Sampling in Three Dimensions: SF^-He 

The SFq He„ system with n = 1 was used to test the rotational sampling 
algorithm. In particular, since for a spherical top we know the exact three- 
dimensional rotational Green's function (Eq. (33) and Appendix A), we can 
check the accuracy of the short-time decomposition of this into three one- 
dimensional rotational factors, Eq. (32). The calculations are carried out in 
the Fixed Laboratory Frame as described in section 4.3, and the SFg-He in- 
teraction is modeled using the anisotropic potential of Pack et al[52]. We 
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Table 2 

Unbiased DMC energies in K for the SFg-He system. The numbers in parentheses 
give one standard deviation statistical errors in the last digits. *The FBR/DVR 
value is converged, i.e., it is stable when the unsymmetrized basis set size is increased 



from 19890 to 33930. 



FBR/DVR 


DMC 


DMC (no rotations) 


-37.14* 


-37.4(2) 


-38.4(3) 



implemented both sampling schemes for this system, i.e., one rotation about 
a randomly chosen axis, corresponding to Eq. (33), and three rotations about 
the space fixed X,Y, and Z axes, corresponding to Eq. (32). Figure 1 (dot- 
ted lines) shows the evolution of the distribution of Euler angles (orientation) 
obtained from the direct sampling of the three-dimensional rotational Green's 
function. It also shows (solid lines) the evolution of the Euler angle distribu- 
tions obtained with the factorization into three independent rotations. Both of 
these methods yield similar distributions that exhibit the expected Gaussian 
form (see Appendix A). The computed ground state energies are summarized 
in Table 2. The value obtained with rotation incorporated is in good agree- 
ment with result obtained using basis set expansion methods. These reference 
calculations employed the Finite Basis Representation - Discrete Variable Rep- 
resentation (FBR-DVR) collocation scheme of Leforestier using a Wigner ba- 
sis set [53]. We also show the energy obtained when the molecular rotation 
is omitted. This is ~ 1 K lower, and is in good agreement with previous 
DMC studies made without rotation [41]. Since the energetic effect of rotation 
is very small for n = 1, i.e., the rotational energy increment in the SFg-He 
ground state is only ~ 1 K, we prefer to investigate the accuracy of the two 
different short time Green's function approximations using an artificial SFg 
having rotational constant Bq = 0.91 cm^^, i.e., 10 times the true gas phase 
value. This increases the variation of the energy due to rotation and amplifies 
the effect of any approximations. Rotating by performing a single rotation 
about a randomly chosen axis yields a ground state energy (Eirot = -33.9(6) 
K) that is consistent with both the value obtained by implementing the three 
rotations about fixed axes scheme (Esrot = -33.8(2) K) and a FBR-DVR refer- 
ence calculation (EpBR = -33.6 K). In contrast, incorrect sampling (omitting 
the a — 3 factor discussed in Appendix A) yields erroneous energies which are 
too low (E = -35.8(7)). This omission is operationally equivalent to reducing 
the rotational constant by a factor of 3. 

5.2 Importance sampling to overcome artificial dissociation in medium-sized 
weakly bound clusters 

The necessity of importance sampling appears already in ground state studies 
of moderately sized, weakly bound clusters, e.g., SFg He„ and HCN He„ with 
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n > 12. Both of these systems require importance samphng to avoid the 
unphysical dissociation of hchum atoms which can result from the occasional 
diffusion of He out to large distances from the cluster. In unbiased DMC, once 



Fig. 1. Euler angle distributions for SFe calculated using the three rotations about 

fixed cartesian axes scheme (solid lines), and using the single rotation about a 
randomly chosen axis scheme (dotted lines). The initial configurations all started 
with the same arbitrary initial molecular orientation. The distributions shown here 
are plotted at three increasing time slices to demonstrate the gaussian spreading. 
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an atom diffuses very far from the weakly bound cluster it has only a very 
small probability to return diffusively to the region of stronger binding, while 
the energetic penalty for motion further away from the molecule becomes 
neghgible. Moreover the configuration space for a motion toward the molecule 
is much smaller than the one for a motion away. This is therefore a metastable 
situation from which it is very difficult to return to configurations of strong 
binding. The result is that helium atoms become "lost" outside the boundaries 
of the simulation box, leading to false dissociation as well as to artificially high 
energies. 

One of the attractive features of unbiased DMC is that it provides the full- 
dimensional wave function. This may be fit to obtain trial functions for im- 
portance sampling, as has been done by a number of authors [32,54,55]. For 
the molecule-helium clusters with n = 1, the helium wave function may be 
squared to obtain the density. However, for n > 1, three-dimensional represen- 
tation of the n-helium wave function is not possible. In order to nevertheless 
provide a visual representation of the wave function, we shall extract an ef- 
fective one-particle wave function by projection, binning independently the 
position of the n heliums in the molecular frame. 

For the spherical top SFg molecule, these calculations are performed in the 
space-fixed, laboratory frame (SF), described in the previous section. For the 
linear HCN molecule, calculations are performed in the Mixed Frame scheme 
of Section 4.1. Rotation around the z-a.xis of the PAF is not allowed {B^ = 0), 
and consequently, the rotation formalism presented there is reduced to two 
angles. 

Simple trial wave functions can be used to avoid artificial dissociation. For one 
rigid molecule, a form of trial wave function possessing the correct permutation 
symmetry is given by the product of pair correlation terms[22]: 

n n 

^TiR, ^) = ll "^He-xiRp, Rx, ^X) n "^He-HeiRp, Rg)- (69) 
p=l pi^q 



In the examples presented here for N < 20, the helium-helium correlation is 
relatively unimportant compared to the helium-molecule term, so for simplic- 
ity we have used ^ He-He{Rp,Rq) — 1- For SFgHei and SF6He2o, the helium 
atoms lie within the first solvation shell. As long as we are interested primar- 
ily in the binding energies here, and not concerned with detailed structural 
information in the angular degrees of freedom, we can adequately describe the 
He-SFe correlations using simple radial wavef unctions, of the form 



"^He-xiRp, Rx, ^x) = exp 



-cr, 



pX 



ar 



(70) 
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where Vpx = \Rp — Rx\, and X=SF6. For this demonstration with n — 20, 
we used a = 0.80 and c = 42000 in atomic units. These parameters were 
obtained from fitting results of unbiased simulations for similar n values. Care 
should be taken not to choose a trial wavefunction that is too restrictive (i.e. 
negligible amplitude in regions that may actually be important). Such choices 
of wavefunction can give rise to instabilities and to erroneous energies. If the 
radial trial function is too narrow, falling off too sharply in the repulsive 
regime close to the impurity, the convergence of the DMC scheme becomes 
more difficult and the energy can become stuck below the true value. We 
found this to be the case when a trial function fit to unbiased wave functions 
obtained for small n {n = 1 with a = 1.5, c = 6000 a.u.) was employed for 
importance sampling of a cluster with significantly larger n {n = 20). This 
emphasizes the importance of choosing a trial function which does not impose 
excessive restrictions on the sampling. 

For HCN, we used the radial part obtained from the high quality anisotropic 
He-HCN trial function derived for use in the fixed node calculations of Section 
V.C. This is of the form 

*t(-Rp, Rx, ^x) = exp cior^^ + C30 In(rpx) - exp(c4o - csorpx) , (71) 



where Vpx = \Rp — Rx\ with X=HCN, and parameters presented in Table 3. 
Table 3 

Parameters (in au) used for the definition of the HCN-He trial wavefunction. 











C2i 


C3i 


C4i 


C5i 


£ = 


0,i = 1 


1.00 


20.07 


-0.051 


-9.18 


5.38 


0.61 


i = 


l,i = l 


1.00 


25.96 


-0.055 


-11.39 


5.66 


0.58 


i = 


2,z = 1 


2.90 


20.17 


-0.031 


-10.14 


6.52 


0.74 


i = 


2,i = 2 


-0.43 


16.38 


0.035 


-9.33 


6.30 


0.85 


i = 


3,z = 1 


4.61 


26.77 


-0.029 


-13.40 


6.68 


0.72 


i = 


3,i = 2 


-6.40 


37.37 


-0.059 


-17.00 


5.44 


0.53 



The effective one-particle helium wave functions obtained by unbiased RB- 
DMC are shown in Figure 2a for SFeHe2o, and in Figure 3a for HCN He2o- 
The "wings" at the outer edge of the binning domains are the signature of 
dissociating helium atoms. They appear here because we put all walkers hav- 
ing a position larger than the domain definition into the last bin at the edge 
of the domain. The average number of artificially dissociated atoms defined 
as those lying at a distance greater than 20 a.u. for both systems, is n^iss ~ 5 
for HCN, and Udiss ~ 1 for SFg. 
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Fig. 2. (a) Radial wavefunction profile for SFg He2o obtained using unbiased DMC 
and projecting the helium wave function into the molecular frame. This is arbitrarily 
normalized as if it were a density, (b) The mixed helium density < ^'t|p|^ > for 
SFg He2o obtained using a radial trial function (see text). 



This dissociation phenomenon also has a pronounced effect on the energy. 
For clusters containing smaller numbers of helium atoms, unbiased RBDMC 
and importance sampled DMC lead to the same energies, as they should if 
both calculations are converged. For n = 20, convergence of the importance 
sampled DMC is straightforward, yielding a ground state energy of -237(4) 
K for HCN and -607 (11) K for SFg. However with unbiased DMC, in sharp 
contrast to this situation, full convergence of the energy is much more difficult 
- even impossible - for this number of helium atoms, because of the artificial 
dissociation of helium walkers, for both HCN and SFg. We can only obtain 
an upper estimate of the ground state energy as a result. For SFe this esti- 
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Fig. 3. (a) One-particle effective ground state wave function for HCN '^He2o, obtained 
using unbiased DMC and projecting the helium wave function into the molecular 
frame. The HCN molecule lies on the z axis, with its center of mass at the origin 
and the H atom on the z > side. The second axis (r) is the distance to that axis. 
The non-physical peaks at the edges of the box come from the binning procedure 
and show artificial dissociation of some helium atoms - see text, (b) Representation 
of the mixed density < > obtained using a radial trial function (see text). 



mate is ~ -564 K, and for HCN it is ~ -196 K. In both cases the unbiased 
DMC value is considerably higher than the fully converged importance sam- 
pled value, reflecting the smaller number of truly bound helium atoms in the 
unbiased calculations. It is significant that this phenomenon is seen not only 
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for the relatively weakly bound HCN molecule in ^He„, but also for the much 
more strongly bound SFg molecule. Our analysis implies that this is a general 
phenomenon which will occur for any species beyond a minimum size which 
appears to be somewhat less than one solvation shell. 

Figures 2b and 3b show the mixed densities < ^r|p|^ > derived from the 
simple radial trial functions with importance sampled RBDMC, for SFg and 
for HCN, respectively. It is evident by comparison with Figures 2a and 3a that 
the use of a simple radial function which provides binding of the helium atoms 
to the embedded molecule very effectively removes the artificial dissociation 
problem, and allows us to fully converge a ground state energy calculation. 
Although we do not evaluate structural quantities here other than this mixed 
density, we note that the use of importance sampling does allow a relatively 
straightforward access to computation of structural quantities using second- 
order estimators [22]. 

Table 4 summarizes the energies and the number of artificially dissociated 
hehum atoms for both clusters. 

Table 4 

Comparison between unbiased and biased RBDMC for SFg and HCN with 20 he- 
lium atoms. No error estimates arc shown for the unbiased results, because the 
convergence with this approach is poorly defined (see text). 





Unbiased energy in K Udiss 


Biased energy in K 


SFq He2o 


-564 ~1 


-607 (11) 


HCN He2o 


-196 ~5 


-237 (4) 



5.3 Excited States 

Our third application addresses the computation of excited states for the He- 
HCN complex. Experiments in doped helium clusters show that the rotational 
spectrum of small molecules inside clusters of ^He possesses the same sym- 
metry as the corresponding spectra in the gas phase, but that the effective 
rotational constants are reduced [56,57]. For HCN this reduction is about 
20%, [58,59] considerably less than the ~ 80% reduction seen for the more 
strongly bound species such as SFg [60]. Table 1 shows that HCN is much 
more weakly bound to helium than SFg, but that in both cases the ground 
state energy of the n — 1 complex lies above the saddle point of the interac- 
tion potential. Analysis of such rotationally excited states within a fixed node 
approximation requires some estimate of the relevant nodal structure. Insight 
into this can be obtained from level assignments made according to simple 
models. The earlier study of SFg in ^Hejv from our group calculated rotational 
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levels corresponding to J ^ j. where J is the total angular momentum, and j 
the molecular angular momentum [17]. For He-HCN, we evaluate here instead 
the excited state which has been assigned by Atkins and Hutson [61], and 
by Drucker et al [62], to the level |j = 1, £ = 1, J = >. In this assignment 
scheme, J is the total rotational quantum number, j the quantum number of 
HCN, which is to a good approximation conserved in the weakly bound com- 
plex, and £ an "orbital" quantum number associated with the rotation of the 
helium around the HCN. For this excited state the fixed nodal approximation 
can be imposed in the molecular frame, as we describe below. Nevertheless, we 
employ here the Mixed Frame implementation of the rotational importance 
sampling algorithm, rather than the Fixed Molecular Frame version. For a 
linear molecule, while we must always perform the rotations in the molecular 
frame, nodal structures may be imposed in either the MF or the SF frame, 
leading to the natural choice of the Mixed Frame implementation as the most 
general in this case. 

In order to obtain a trial wave function for the fixed node calculations, we 
first performed pure DMC runs which allow us to compute the ground state 
eigenfunction (see Figure 4). This was fit to an analytic form ^'T(r, 6*), where 



Fig. 4. Ground state helium wave function obtained by pure DMC propagation pro- 
jected into the molecular frame for He-HCN. The orientation of the HCN molecule 
is the same as in figure 3. 



r and 6 arc the usual Jacobi coordinates of the He with respect to the center 
of mass and orientation of HCN, to obtain a high quality trial function. We 
use the following Legendre polynomial expansion: 
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3 

e=o 

2 

"^iir) = ^cii exp [cur'''^' + c^i ln(r) - exp(c4i - c^ir)] 

i=l 

with {ai,Cki,k = l,...,5;i = 1,2} the parameters of the fit presented in Table 
3. Limiting the sums to £ < 3 and i <2 does not noticeably affect the precision 
of the fit. 

Such a high quality trial function is necessary for the fixed node calculations of 
the He-HCN complex described here. In order to construct a trial excited state 
function for a fixed node calculation, we combined this high quality ground 
state trial function with a simple function that imposes the nodal constraint. 
This is chosen here by comparison with the exact excited state function, which 
was obtained by repeating the collocation calculations of Aktins and Hut- 
son [61] to compute and analyze the wavef unction for this level. The energy of 
this level is —8.5 K, and a contour plot of the corresponding wave function 
is shown on figure 5. The structure of this function implies a strong internal 

Fig. 5. Excited state |110) wavefunction obtained by DVR calculation represented 
as a function of the Jacobi coordinates r in A and 9 in degrees. The contour lines 
are evenly spaced ranging from -0.5 to 0.4 



bending character to this excitation. The excited state function possesses a 
single nodal surface, which is seen to be approximately r-independent. This 
motivated us to employ a fixed nodal surface defined by cos(^ -|- x), where x is 
a parameter and 6 is the internal Jacobi angle of the cluster i. e. the trial func- 
tion used is thus "^rif^, 9) = ^^(r, 6) cos{6 + x)- Since the calculation is carried 
out in the mixed frame formulation, this trial function is implicitly dependent 
upon the orientation of the rigid body in the space fixed frame, even though 
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the nodal structure is imposed in the molecular frame. Comparison of the en- 
ergy obtained for calculations restricted to each of the two sides of the nodal 
surface allows us to then obtain the optimal value of the parameter x- This is 
arrived at when these two energies are equal. A similar procedure was used by 
McCoy and co-workers in an interactive sense to optimize nodal surfaces [32]. 
For both the ground state and excited state computation, we use an ensemble 
of approximately 5000-6000 walkers. We perform a block averaging in order 
to reduce the correlation between successive steps. After equilibration of the 
ensemble, we perform 1 run of 800 blocks, each of which consists of 150 time 
steps with At = 10 au. We reject all moves corresponding to a crossing of the 
nodal surface, keeping the corresponding walker at its former position [39]. 

The optimal value of x was found to be x = 15.65 degrees, which yielded a 
common energy of -8.4(2) K on both sides of the resulting nodal surface. This 
fixed node value is in excellent agreement with the vahie -8.48797 K (-5.89897 
cm~^) obtained from the collocation method in Ref. [61], especially when the 
approximate nature of the nodal constraint is taken into account. 

When studying clusters with larger numbers of helium atoms using unbiased 
RBDMC, we noticed that the 9 dependence apparent in Figure 4 tends to 
smooth out [63]. This suggests that for larger clusters, it might be adequate 
for some purposes to reduce the trial function to the radial term only, e.g., for 
studying ground state energies (Section V.B). 



6 Conclusions 

We have derived an algorithm for Rigid Body Monte Carlo with importance 
sampling for all degrees of freedom, including both translation and rotation, 
and provided a rigorous derivation of the short time rotational Green's func- 
tion and its associated quantum forces. Three possible implementations of the 
algorithm were presented, which differ according to the combination of frames 
used for rotational and translational moves during the propagation. The most 
general implementation scheme is the mixed frame implementation in which 
translational moves are made in the laboratory reference frame, while rota- 
tional moves of all rigid bodies are made in the principal axis frame of each 
monomer. This allows full importance sampling calculations of molecular clus- 
ters to be performed. 

We have presented several applications of this importance sampled Rigid Body 
Diffusion Monte Carlo. First, we demonstrated the correctness of the short- 
time factorization of the rotational Green's function by comparison of the 
resulting product of one- dimensional rotations with the well known three- 
dimensional formulation for a spherical top. Second, we showed that impor- 
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tance sampling of translational degrees of freedom is essential to avoid non- 
physical dissociation of weakly bound helium atoms in doped quantum clusters 
of helium when more than a few helium atoms are present. Third, we showed 
how excited states can be easily accessed with this algorithm using the fixed 
node approximation and a trial wavefunction with an implicit dependence 
on the orientation of the rigid body. The same algorithm may now be ap- 
plied to study excitations beyond the fixed node approximation, e.g., using 
the POITSE approach[31]. 
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8 Appendix 



8.1 Sampling of the Spherical Top Green's Function 



In order to provide a sampling of the operator 



(72) 



where d= 1/27, according to the configuration representation of Eq. (33), one 
must choose both an arbitrary axis of rotation, and an angle of rotation. In 
order to implement this, it is necessary to first specify the frame of reference of 
the operator L, i.e., its axis of quantization for L^. This is evident on expansion 
of the operator in the angle space. 



where 1^ is a unit vector in an arbitrary direction, and a is a constant to be 
determined. Expanding both sides of Eq. (73) to order (Ar)^, we find 

1 - dL'^Ar = 1 - ^ y dn{L ■ Clf. (74) 
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The integral can be evaluated once the direction of quantization of L is estab- 
lished. Choosing this for simplicity as z, so that L-Cl = Lcos6, and evaluating 
the integral over the spherical angular coordinates dfl = d4>d9, leads to the 
equality 

1 - dL^Ar = 1 - d-L^Ar, (75) 
3 

from which we conclude that a — 3. Sampling of Eq. (73) can then proceed 
by a) choosing a random vector Cl, and then b) sampling the angle of rotation 

by choosing a value from a gaussian of width y^3 x (2d At). This procedure 
provides access to the full angular space in Eq. (73). When we are dealing with 
small rotations, this approach gives rise to similar orientational distributions 
to the scheme in which 3 small rotations are performed about the 3 cartesian 
axes (see Figures 1). This similarity will continue to hold as long as the rotation 
C is small enough such that the single rotation about the arbitrary axis can 
be approximately factored into a product of three cartesian rotations, i.e., 

exp (iC J • n^) ~ exp (i^Jx sin 9 cos 0) exp [iCJy sin 9 sin (p) exp {iCJz cos 9) . (76) 

Then the corresponding one-dimensional rotation around, e.g., the x axis, is 
made with the small angle C sin 9 cos (p. 



8.2 Commutation Relations 



Expanding the standard commutator 



d 



dx 

we arrive at 



df{x) 
dx 



(77) 



dx dx dx 
h "'•'^ ' dx 



(78) 
(79) 



Generalizing the one-dimensional function f{x) to the three-dimensional F{R) - 
fxi + fyi + fzk-i and evaluating Eq. (79) for each degree of freedom, leads then 
to the desired relation 



F- V = -P-F- V-F. 

h 



(80) 
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8.3 Computational Scheme 



Start with an initial ensemble of walkers (e.g. distributed according to the trial 
wavefunction). Compute the local energy (also used as the starting reference 
energy) for the initial ensemble and then propagate the ensemble in imaginary 
time by repeating the following steps for each walker (configuration). 

(1) Move each host and impurity atom according to equation (60). This move 
involves a Gaussian random number and a quantum force. 

(2) Each move is then accepted or rejected with acceptance probability given 
by mm(l, H^((5', Q)) where 



|^T(QOpg(Q^^Q,AT) 



This involves having to recompute the quantum force after the particle 
has been moved. 

(3) Rotate the impurity atom according to equation (61) 

(4) Again, this move is accepted or rejected using the above criterion. 

(5) Determine effective time step, Are//. See Ref. [39]. One has to average 
the effective time step over the various kinds of moves since in general, 
each kind of move has a different diffusion constant. 

(6) Compute the new local energy, El{Q') (equation (37)) 

(7) Compute branching factor, M: 



M = int exp 



+ A 



where A is a random number that is uniformly distributed over (0,1). 
If M = 0, this walker is destroyed, otherwise M — 1 duplicates of the 
configuration are added to the ensemble. 
(8) To maintain a stable population size (N(r)), update the reference energy 
according to 

where the population control parameter, a is chosen to be as small as 
possible (to avoid biasing the results) but large enough that the popula- 
tion size stays within acceptable hmits. N* is the desired population size 
and is often chosen to be the starting population A'"(0) or N{t). 

Of course, many variants of the above procedure will also work. In the HCN 
studies, for example, all of the translation moves and rotations are performed 
simultaneously and then accepted or rejected as a single move. When studying 
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excited states using the fixed-node approximation, one also adds the additional 
constraint that moves crossing the nodal surface are automatically rejected. 

The statistical error is estimated using block averaging, where the propagation 
is split into Nj, blocks of Ng steps where A^, is longer than the correlation 
length [47]. In the literature, various methods are used to compute averages 
and errors from this blocked data: 

(1) Take 1 data point per block, and compute the average and standard error 
from these Nb data points 

(2) Same as above but use the block average values for the A*";, data points 

(3) Use all the data points in computing the mean. The standard error for 
each block, a^, is computed from the Ng values within each block. These 
are averaged to determine the reported error: 

In our examples, we report the largest of the above estimates (usually close in 
value) . 
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